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Abstract 

We describe the implementation of an angular power spectrum estimator in the flat sky approximation. POKER (P. 
Of k EstimatoR) is based on the MASTER algorithm developped by Hivon and collaborators in the context of CMB 
anisotropy. It works entirely in discrete space and can be applied to arbitrary high angular resolution maps. It is 
therefore particularly suitable for current and future infrared to sub-mm observations of diffuse emission, whether 
Galactic or cosmological. 
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1. Introduction 

Whether it is due to Galactic dust or synchrotron, to 
cosmological backgrounds such as the Cosmic Microwave 
Background (CMB) or to the Cosmic Infrared Background 
(CIB), that traces the integrated radiation of unresolved 
galaxies, diffuse emission is omnipresent in infrared and 
millimetric observations. The angular power spectrum of 
this radiation is one of the main tools used to constrain 
the structure of the interstellar medium, the clustering of 
IR galaxies (CIB), or the cosmological parameters (CMB). 
In short, its estimation requires to Fourier transform the 
image and to average the modulus square of the Fourier 
amplitudes into frequency bins. However, the image has 
both boundaries and often masked regions (e.g. to re- 
move bright point sources) that induce power aliasing 
and biases the estimation of the power spectrum if not 
accounted for properly. The effect becomes quite signifi- 
cant when the signal has a steep power spectrum such as 
similar to that measured for Galactic cirrus emission 



(jMiville-Deschenes et al.l 120071) or even steeper than k 4 
as for CMB ani sotropy on angular scales smaller than a 
few arcmin (e.g. iReichardt etall 120091: iBrown et ai]l2009lh 
To account for non-periodic boundaries. iDas et al. r (|2009h 
proposed an original apodizing technique that helps us to 
deconvolve the estimated power spectrum from that of the 
observed patch boundaries. These authors also mitigate the 
impact of holes by applying a pre- whitening technique to 
data in real space. 

In the context of CMB anisotropy, iHivon et ail ([20021) 
developed the MASTER method that allows us to correct for 
mask effects on the output binned power spectrum. They 
analyze data accross the full sky and account for the sky 
curvature. Instead of classical Fourier analysis, they project 
the data onto spherical harmonics and go through the al- 
gebra of pseudo angular power spectra (see Sect. 0). This 
idea has been successfully used in several experiments (e.g. 



Ide Bernardis et al.l 120001: IBenoit et all 120031) and is also 
the basis of more re fined algorithm s used in e.g. WMAP 
(jHinshaw et al.ll2003l ) and Archeops dTristram et al.ll2005l ). 
However, direct use of MASTER in the context of infrared 
observations with a resolution of typically a few arcsec re- 
quires us to estimate Legendre polynomials up to orders 
£ of 10,000 or more for which current recurrences and in- 
tegration methods are numerically unstable. Other tech- 
niques developed in the context of CMB anisotropy such as 
maximum likelihood estimation dBond et al.l 119981 ) could 
be transposed to high angular resolution maps, but the nu- 



merical cost 



is prohibitive for common applications 
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when the analysis pipeline requires Monte Carlo simula- 
tions. 

This paper aims to transpose the pseudo-spectrum ap- 
proach pioneered by MASTER to high angular resolution 
observations and a classical Fourier analysis in the con- 
text of the flat sky approximation. Its originality compared 
to other approaches is that it works exclusively in discrete 
space and therefore avoids the complexity of resampling the 
data and integrating Bcsscl functions. Our algorithm was 
nicknamed POKER, for "P. Of k EstimatoR". The paper 
is organized as follows. Section [2] provides the definitions 
and algebra involved in POKER and Sect. [3] shows its appli- 
cations to simulations of various astrophysical components 
spectra and a complex mask. Detailed derivations of our re- 
sults arc presented in the appendices. Although this work 
focuses on temperature power spectrum estimation, we also 
show in Appendices IBI and [Clhow the formalism can be gen- 
eralized to the case of polarization. 

2. Power spectrum estimation for an incomplete 
observation of the sky 

We first briefly state the limits of the flat sky approxima- 
tion, then recall the definitions of both the power spectrum 
and pseudo-power spectrum of data in the context of con- 
tinuous Fourier transforms. Finally, we consider their coun- 
terpart in discrete space and the implementation of POKER. 
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2.1. Flat sky approximation 

Projecting an observed fraction of the sky onto a plane 
rather than a sphere alters the image properties in a 
way that depends on the specific reprojection scheme (e.g. 
gnomonic, tangential, cylindrical etc). The angular power 
spectrum of the data as measured on the projection plane 
therefore differs from that on the sphere. Two comments 
can be made at this stage. First, as long as the observation 
patch does not span more than a few degrees, as in most 
infrared experiments, the distortions are very small. For a 
gnomonic projection for instance, an observation point at 
an angular distance 9 from the map center (the point tan- 
gent to the sphere) is projected at a distance tan(#) rather 
than 9. The relative difference between the two is only 
2.6 10~ 3 for a map of 10 degree diameter i.e. the projected 
map is stretched by 46 arcsec in each direction. As long 
as this remains small compared to the angular scales over 
which the angular power s pectrum is estimat ed, the distor- 
tion can be neglected (e.g. lPrvke et alj |2009). In this work, 
we stay well above this limit by considering square maps 
of only 5x5 deg 2 and 2 arcmin resolution and therefore 
neglect distortion effects. Second, the impact of the projec- 
tion on the estimation of the power spectrum is equivalent 
to a transfer function that can be estimated and corrected 
for when dealing with the convolution kernel (Eq. l2"Tj) , as 
detailed in Sect. 12.51 Both of these reasons allow us to es- 
timate the angular power spectrum of a map built in the 
flat sky approximation limit. 



size of the observation patch. This is accounted for by a 
weight function W r applied to the data. Its most simple 
form is unitary on the data, zero outside the observation 
range or where strong sources are masked out. More subtle 
choices such as inverse noise variance weighting or apodiza- 
tion (cf. Sect. [3]) are usually used. Instead of the true Fourier 
amplitudes, we are then bound to measure the amplitudes 
of the masked data, a.k.a. the pseudo-amplitudes 

f fe = / drT r W r e- lkr . (6) 
Jr 2 

Equation (0| applied to the pseudo-amplitudes gives the 
lD-pseudo-power spectrum 

/>oo 

P fe = / kxdkx K kkl P kl , (7) 
Jo 

where K kkl is the mixing matrix that depends on the 
weighting function W r . To determine the signal power 
spectrum, we need to solve this equation for P kl . A de- 
tailed derivation o f the analytic solution can be found in 
iHivon et~ al. (2002). The impact of the instrumental beam, 
the pixel window function, the projection algorithm, the 
scanning, and the data processing can be accounted for in 
this formalism as transfer functions and incorporated in the 
definition of K kkl (Sect. |2~5"1) . In the next two subsections, 
we focus on mask effects and the global picture of our al- 
gorithm. 



2.2. Continuous Fourier analysis and masked data 

On a flat two-dimensional (2D) surface, a scalar field T r 
depending on the direction of observation r is represented 
in Fourier space by 

T k = [ drT r e~ ihr , (1) 
Jr 2 

Tr = f j^T k e* kr . (2) 
Jr 2 (2tt) 2 

For a random isotropic process, the 2D power spectrum V k 
is defined as 



2.3. Discrete Fourier analysis and POKER 

Any data set is by construction discretely sampled. 
Computing the quantities defined in the previous section 
requires mathematical interpolation and/or resampling of 
these data and appropriate integration tools, especially if 
the underlying data power spectru m is steep as for Galactic 
cirru s for which P(k) cx fc -3 (iMiville-Deschenes et al.l 
120071 ). Rather than dealing with these difficulties, we keep 
the native pixelized description of the data and work com- 
pletely in discrete space. We use the discrete fourier trans- 
form (hereafter DFT) as provided by data analysis software. 
For a map of scalar quantity D^ v and size N x x N y pixels, 
it is defined as 



(T k T*,)=V k S k ^, (3) 

where the brackets denote the statistical average. We de- 
note by one-dimensional (ID) power spectrum its azimuthal 
average 

P k = ^ jj d6T k T* h , (4) 

where P k is the physical quantity of interest which we want 
to reconstruct. It is the Fourier transform of the two-point 
correlation function. If the process is isotropic, the 2D and 
ID-power spectra are related by 

(T k n,) = V k 5 k _y = {2TrfP k 8 k _ k ,, (5) 

In the following, we neglect the ID or 2D qualifiers to 
improve readability and in both ID and 2D cases use the 
term power spectrum unless the difference needs to be em- 
phasized. In practice however, the integrals of Eqs. (fTJ [5]) 
cannot run up to infinity simply because of the limited 



mn ~ N X N V ^ ' [ ' 

Dp,v = £) mne +2iir{ f j,m/N I +vn/N y ) ^ fg^ 

m,n 

Throughout this work, although we denote quantities in 
direct and Fourier space by the same name, Greek indices 
denote pixel indices in real space whereas roman indices 
refer to amplitudes in Fourier space. Unless stated other- 
wise, sums over /x and m (rcsp. v and n) run from to 
N x — 1 (rcsp. N y — 1), A9 is the angular resolution of the 
map in radians. For a given wave- vector k mn , labeled by 
the m and n indices, its corresponding norm is denoted by 
k mn = (2?:/A9)^(m'/N x ) 2 + (n'/N y ) 2 with m! = m (resp. 
n') if m < N x /2 and m' = N x - m if m > N x /2. This con- 
vention ensures that on small angular scales k matches the 
multipolc £ used in the description of CMB anisotropy. The 
Nyquist mode is tt/ A9. 
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Figure 1. Map of the Fourier modes of the worked exam- 
ples of Sect [3] The inner circle delimits the Nyquist range. 
Modes that lie on the outer circle are examples of modes of 
larger modulus than k^yquist- For these modes, not all di- 
rections are sampled in the Fourier plane (dashes represent 
the missing modes). 



It is well known that the DFT slightly differs from the 
theoretical continuous Fourier transform, hence D mn does 
not strictly equal Xfe m „ . In particular, the DFT deals with 
amplitudes for modes k mn larger than the Nyquist mode 
ir/A9 and in some directions for 9 mn only (see Fig.[T]). It is 
therefore not possible to integrate Eq. Q on the full range 
9 G [0, 2tt] for such modes and so, the ID power spectrum 
is undefined outside the Nyquist range. In the following, we 
therefore restrict ourselves to the Nyquist range for power 
spectrum estimation. We note, however, that mathematical 
sums implied in the following may still run over the full 
range of pixels or DFT amplitude indices. 

The direct DFT of the masked data results from the 
convolution of the DFT amplitudes by a kernel K^f 1 , that 
depends only on the mask DFT amplitudes (Appendix [5J . 
If the data D consist of signal T and noise TV, we have 



(|Ann| 2 ) 



E 

m'n' 



\T„ 



+ (N rnn ) , 



(10) 



which is the transcription in discrete space of Eq. (0 . 

The rapid oscillations of the convolution kernel intro- 
duce strong correlations between spatial frequencies and 
make its inversion numerically intractable. (PseMC?o-)Power 
spectra are therefore estimated on some frequency band- 
powers (labeled b hereafter). The binning operator reads 



*w if i.b 
Eh 11 K low 



^ kmn ^ & 



,6+1 



otherwise 



(11) 



where kf ow is the mode of lowest modulus that belongs to 
bin b and St is the number of wave vectors k mn that fall into 
this bin. The reciprocal operator that relates the theoretical 
value of the ID binned power spectrum Pj to its value at 

km/n IS 



Q 



b 

mn 



if k b 
11 How 



_^ ">mn ^ ™\ 0I 

otherwise 



6+1 



(12) 



Although not strictly required, results may be improved 
when the spectral index /3 is chosen so that k^Pk is as flat 
as possiblqE In the case of the cosmic infrared backg round 
anisotropy, (3 ~ 1 (jPlanck Collaboration et all 1201 it ). The 
binned pseudo-power spectrum is therefore given by 



Ph 



E * 



b I T mn I 



(13) 



and the data power spectrum is related to its binned value 
Pb via 



IT / A 2 

-L m'n' 



Qrn'n'Pb' ■ 



With such binned quantities, Eq. (|10[) reads 
(? b ) - E M w p v + (N b ), 

V 

where 



M w = 



E E 



K b \ K m,m' 



2 O b 



(14) 
(15) 

(16) 



An unbiased estimate of the binned angular power spec- 
trum of the signal is thus given by 



n-E^'fA'-^ 



(17) 



It can indeed be easily verified that (Pb) = Pb- Uncertainties 
in Pb come from sampling and noise variance that are esti- 
mated via Monte Carlo simulations as described in the next 
section. 

2.4. Statistical uncertainties 

Statistical uncertainties in Pb come from signal sampling 
variance and noise variance. They are estimated via Monte 
Carlo simulations. For each realization, a map of sig- 
nal+noisc is produced (see Sect. I2.6[) and treated in the 
same way as described for the data in the previous section 
to give an estimate, P\, with the same statistical proper- 
ties as that of the true data. Altogether, these simulations 
provide the uncertainties in our estimate. The covariance 
matrix of Pb is 



C bb ' = ({Pb- (PbUc) [Pb' - (Pb')Mc) ) 



'MC 



(18) 



where (-)mc denotes the Monte- Carlo averaging. The error 
bar in each Pb is 

<J Pb = VCb~ b , (19) 

and the bin-bin correlation matrix is given by its standard 
definition 



-bb' = 



b'b 1 



(20) 



1 In the case of CMB, /3 ~ 2 is the equivalent of the standard 
{I + 1) prefactor that flattens the spectrum up to I ~ 2000. 
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Figure 2. Top: Mask applied to the simulated data. This 
mask is 1 where data are available and outside the obser- 
vation patch and where bright sources have been masked 
out. Bottom: Same mask with apodizcd boundaries but still 
with the same number of masked pixels. 



2.5. Beam, map making and transfer functions 

The above sections describe the main features of the al- 
gorithm that provides an unbiased estimate of the power 
spectrum of data projected onto a map. This power spec- 
trum may however differ from the signal power spectrum. 
The map making process may indeed alter the statistical 
properties of the signal, together with data filtering and 
the convolution by the instrumental beam. For instance, 
the pixelization caused by the map making is equivalent to 
a convolution in real space by a square kernel and there- 
fore translates into a multiplication in Fourier space by a 
factor sin c . In the case of CMB anisotropy for which power 
spectrum estimation has been most extensively studied, to 
a good approximation, the transfer function of the map 
making and the data processing together reduces to a func- 
tion of k that multiplies the data power spectrum P^. The 
determination of Fk is performed by a set of Monte Carlo 
simulations of data processing and map making. The beam 
smearing effect is also described by a multiplicative func- 
tion Bk- In the present framework, it is possible to be even 
more precise and to account for the exact beam shape and 
orientation because the beam can be completely described 
by its Fourier coefficients B mn rather than its approximated 
annular average Bk ■ This may be of particular relevance for 



small fields over which the scanning directions arc approxi- 
mately constant and increase the effect of beam asymmetry. 
The map making together with the filtering transfer func- 
tion is also likely to be more accurately represented by a 
function of both Fourier indices F mn , such that Eq. (fTU)l is 
given by 

(\D mn \ 2 ) = \K^U 2 F mlnl B m , nl P km , n , 

m'n' 

+{N mn ). (21) 

These new contributions can be incorporated in the defi- 
nition of the convolution kernel K^^i such that no addi- 
tional modification of the algorithm is needed from Eq. ([TUll 
onward. 

2.6. Algorithm 

An outline of the algorithm is presented in Fig. [5] To per- 
form the complete process of power spectrum and statistical 
uncertainty estimation, one needs: 

(a) A tool to simulate the sky that is observed by the in- 
strument given an input angular power spectrum. 

(b) A tool to simulate the instrument observations given 
the simulated sky. 

(c) The data processing pipeline that derives from these ob- 
servations the map from which the user can then esti- 
mate the angular power spectrum. The pipeline includes 
optical beam smearing, time domain filtering, data flag- 
ging, map making etc. This tool is required to determine 
the transfer function of the map making and data pro- 
cessing (Eq. l2"Tjl . 

(d) A tool to compute the power spectrum of a 2D map and 
to bin it into a set of predefined bins with a weight that 
may be a function of k. 

(e) A tool to compute M^w, which involves the computa- 
tion of the convolution kernel K^ 1 ^,. This is actually 

the longest part of the process because it is a N^ ix op- 
eration, but it only needs to be done once. 

All these tools are provided in the POKER librarjffl. The 
algorithm can be summarized as follows: 

1 . Insert the observed sky patch of size N x x N y pixels into 
a "large patch" (N' x x N' y ) and padd it with zeros. This 
will allow for the correction of aliasing by scales larger 
than the observed sky. The size of the patch and the 
zero padding that should be used have both to be de- 
termined by the user. A factor from 1.2 to 2 is enough in 
most cases. A compromise must be chosen between the 
uncertainty on large scales that the user tries to esti- 
mate and the uncertainty associated with the unknown 

2 http://www.ias.u-psud.fr/poker This libr ary m akes use 
of some of the HEALPix programs (|Gorski et alJfl99Sh . 

3 The sky simulator provided in the library works directly in 
fiat space. Indeed, we cannot anticipate neither on the user's 
favorite reprojection scheme from the sphere to the plane nor on 
his/her data processing pipeline and map making. Furthermore, 
simulating directly a fiat sky of a known power spectrum is a 
convenient first stage to optimize the binning and apodization 
scheme before running the full fledge Monte Carlo simulations. 
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power in these large scales that needs to be assumed for 
the simulations. It is also possible to apodize the obser- 
vation patch to limit large-scale aliasing (see Sect. |3] for 
more details) and improve the bin to bin decorrelation 
at high k. 

Define a binning for the estimated power spectrum 
on the large patch. Typically, modes sampled by the 
data set are the DC level and modes between k m i n — 
2tt / A8 / ma.x(N x , N y ) and the Nyquist mode k c = 
n/A9. The minimum bandwidth of the bins may be 
chosen as ~ 2fc„ UIl . 

Determine the noise pseudo-power spectrum - (Nb) of 
Eq. (fT5j) . If it cannot be determined analytically, per- 
form a set of Monte Carlo realizations of noise-only 
maps (with (a)) and compute the power spectrum with 
(b) of the masked maps inserted into the large patch. 
The average of these Monte Carlo realizations gives 

(N b ). 

Run a set of noise-free simulations of the observed and 
rcprojcctcd sky to determine the transfer function F mn 
of the data processing pipeline - Eq. ([2"T1) . 
Compute M bb > with (c) - Eqs. ([TT1 [121 1TB ]) . This opera- 
tion scales as but it only needs to be done once. The 
implementation proposed in the POKER library can be 
run on a multiprocessor machine. 

Compute the pseudo-power spectrum of the masked 
data on the large patch Pb with (b) - Eq. (jA.ip . 
Apply Eq. (fTT)) to obtain the binned power spectrum of 
the data Pb- The resolution of this equation can be done 
with any suitable method of linear algebra. Note that 
Mbb' can be rather small and its inversion straightfor- 
ward with standard numerical tools such that Eq. (fTT)) 
can be computed as is. At this stage, it may be useful to 
discard the first bin of the matrix, which describes the 
coupling of the DC level of the map to the mask and is 
therefore irrelevant for a power spectrum analysis but 
tends to alter the conditioning of Mbb' ■ 
Determine the statistical error bars associated with this 
estimate. For that, perform a set of Monte Carlo real- 
izations of signal+ noise. The input spectrum required 
for these simulations can be a smooth interpolation of 
the binned power spectrum determined at the previous 
step. For each realization, compute the pseudo-power 
spectrum (using (b) on the masked data embedded in 
the large patch), subtract (Nb), and then solve Eq. (fTT)) . 
This provides a set of random realizations of Pb- The er- 
ror bars and the bin-to-bin covariance matrix are then 
given by Eqs. (|18I19D . 



3. Worked example 

POKER was applied to real data to measure the cosmic 
infrared background anisotropy i n the Planck-HFI data 
(jPlanck Collaboration et ahl 1201 lh . The whole data pro- 
cessing and how it is accounted for with POKER is described 
in detail in that paper. Wc do not replicate this analysis 
here but instead present complementary examples on sim- 
ulated data with steeper power spectra and a more complex 
mask. It is in this context that mask aliasing effects are the 
strongest. 

We assume that the observation patch is a square of 100 
pixels of 2 arcmin side. These parameters arc chosen so that 




Lon (deq) 




Lon (deg) 



Figure 3. Typical maps of dust with a k~ 3 power spectrum 
(top) and CMB temperature (bottom). The square is the 
outline of the observed patch. It is extracted from a larger 
simulated map to ensure non-periodic boundary conditions. 
Masked data appear in white. 



they sample a range of angular modes over which the CMB 
temperature power spectrum exhibits peaks and a varying 
slope. Note however that the map resolution can be arbi- 
trary high because fast Fourier transform algorithms work 
in dimensionlcss units. To force non-periodic boundary con- 
ditions, we extract the patch from a map that is 50% larger 
and the simulation is performed on the latter. Finally, we 
draw random holes across the observation patch to mimic 
point source masking. We consider two types of signal. In 
the first case, we assume that the data are represented by a 
pure power law spectrum A; -3 typical of Galactic dust emis- 
sion. In the second case, we assume that the data are CMB 
with a standard ACDM power spectrum. At these angular 
scales, the slope of the CMB power spectrum varies from 
~ k~ 2 to even steeper than k~ 6 and exhibits oscillations. 
Fig. [3] shows an example of these simulated data. The re- 
sult of POKER applied to each case is presented in Figs. @] 
and[SJ In the case of dust, we choose a binning index (3 = 3 
as defined in Eq. ([T2]) . and in the case of CMB, we make 
no assumptions, i.e. we choose (3 = 0. We also show what 
the direct Fourier transform of the observed patch without 
any further correction would give to illustrate the magni- 
tude of the effect corrected by POKER. Note that this refer- 
ence estimate labeled "naive P(k) n is not the pseudo-power 
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Figure 4. Dust (k~ 3 ). Comparison between the input the- 
oretical power spectrum (black) and the average result of 
POKER (red) applied to 500 signal+noise simulations. The 
"naive" approach (blue, see text) is also shown for refer- 
ence. Error bars in the top plots are those associated with 
the data (i.e. those of a single realization). The square 
line shows the binned theoretical power spectrum to which 
POKER's average result should be compared. The bottom 
plots shows the ratio of the reconstructed binned power 
spectrum to the input theoretical binned power spectrum 
(the bias) and the displayed error bar is that of the av- 
erage of the Monte Carlo realization (in other words: the 
error bar of the top plot divided by \/500). These plots al- 
together show that POKER is unbiased. The mask used in 
this case is that of the top plot of Fig. [3J with 1 where data 
are available, elsewhere. 



spectrum of the data in the sense of Sect. [2] Indeed, it is 
not computed on the whole map from which the observa- 
tion patch is extracted and padded with zeros. The bottom 
plots of Figs. [4] and [5] show the bin-to-bin correlation ma- 
trix of each estimate. In the case of dust, the correlations 
are small (~ 15%). This is not the case for the simulated 
CMB, for which there is strong bin to bin correlation, al- 
though the power spectrum remains unbiased. This corre- 
lation is due to large-scale aliasing induced by the holes in 



the mask and show up so significantly because at high k, 
the CMB spectrum is very steep. A way to improve on this 
is to apodize the mask around the edges and the holes left 
by point-source masking (Fig. [5]). In this work, we simply 
use a Gaussian kernel with a FWHM of twice the map res- 
olution to smooth the edges. The same analysis as before 
is performed with this mask and results are presented in 
the right hand side of Fig. [SJ On this occasion, the bin-to- 
bin correlation is significantly reduced, albeit the sampling 
variance is slighty increased at low k owing to the effec- 
tive reduction in the observation area. A more efficient way 
of per forming this apodization is described in iGrain et al.l 
(|2009t ). Finally, on larger angular scales, there is a slight 
bias in the recovery of the CMB power spectrum. This 
does not however occur in the case of dust, because for 
a pure power-law spectrum, Eq. (|14[) is then an equality. 
This is no longer the case for a CMB spectrum whose av- 
erage slope varies with k and exhibits peaks. No binning 
could faithfully represent such a spectrum. However, the 
remaining bias is negligible relative to the statistical error 
bar in the data. If we force the simulated CMB to have a 
constant power spectrum over frequency bins, the recovery 
is unbiased. There is no general prescription regarding the 
definition of the binning and the apodization. They must 
however be chosen with care because the bin-to-bin resid- 
ual correlation may lead to residual ringing (mask aliasing) 
in the data power spectrum (considered as a single random 
realization), even if the estimator is, on average, unbiased. 



4. Conclusion 

We have developed a tool that provides an unbiased esti- 
mate of the angular power spectrum of diffuse emission in 
the flat sky approximation limit, for arbitrary high resolu- 
tion and complex masks. POKER corrects for mask alias- 
ing effects, even in the context of steep power spectra and 
provides a way to estimate statistical error bars and bin- 
to-bin correlations. It complements tools developed in the 
cont ext of spherical sk y and potentially full sky surveys 
(e.g. iHivon et al.ll2002t ) but at the moment for lower angu- 
lar resolutions. POKER also comple ments other meth ods in 
the flat sky approximation such as iDas et all (|2009() . 

POKER can readily be generalized to polarization 
power-spectra estimation. To date, exp eriments tha t have 
measured polarized diffuse emission (|Kovac et aT 
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Chiang et al 



2003 
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2002 



2008 



l2010t Bierman et al.l 20 111) have been closely 



related to CMB experiments and studied observation 
patches of a few to a hundred percent of the sky and an- 
gular resolutions larger than a few arcmin. Optimal tools 
have been developed to measure the polar i zation power 
spectra in this context (IChon et al.l l200i ISmithl l200l 
ISmith fc Zaldarriaeal l2007t IGrain et al.l 120091 " and refer- 
ences therein) and it is unlikely that POKER will bring 
something significantly new to the analysis of these obser- 
vations. It is however expected that smaller, deeper, and 
higher-resolution polarized surveys will happen in the fu- 
ture, for which POKER might be an interesting approach. 
One of the main features that should then be addressed is 
the ability of POKER to correct for E — B leakage. Although 
we postpone the detailed studies of POKER's properties re- 
garding polarized power spectra estimation to future work, 
we provide the formalism in Appendices [B] and [C] for the 
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Input Data - 
Input Data + noise - 
Naive P(k) - 
POKER - 




Input Data - 
nput Data + noise - 

n °7o p kS : 



Figure 5. Left: same as Fig. @] in the case of CMB. POKER's estimation is unbiased but mask aliasing induces strong 
correlations between bins at high k where the power spectrum is very steep. Right: This time, the mask with apodized 
boundaries (bottom of Fig. [2]) is used. High k bin-to-bin correlations are significantly reduced. Apodization however 
reduces the effective observed fraction of the sky and therefore slightly increases error bars at low k compared to the 
plots on the left. Note that although apodization also improves the "naive" estimate, it remains not compatible with the 
input spectrum for almost every bin and to more than the la error of a single realization. 



sake of completeness. All the software used in this work is 
publicly available^. 
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38 




I V 




- 




Binning 



DFT, Eqs. 11, 12, 13 

r 



pi p2 p,V 




Figure 6. Schematic flow chart of POKER. 

Appendix A: Mask convolution kernel for temperature 

If not specified, sums run from to N x — 1 and to N y — 1. The pseudo-Fourier coefficients of the weighted data are given by 



m " ~ N X N,- ^ M " 



ST ST T „2iTr(iim 1 /N 1: + vn 1 /Ny) ST m „2in(nm 2 /N x + vn 2 /N y ) -2iir(ltm/N x + vn/N y ) 

2__, 2—1 lm l n l e 2—1 vv m 2 n 2 e e ! 

m 2 ,n 2 

2i7r[^(mi iV x +m 2 — m) /iVa, -\-i/(n 1 +n 2 — n)/N y ] 



N X Ny 



1 



N X N, 



52 52 T ™^ W < 



m 2 n 2 c 



raj ,ni ,m 2 ,n 2 ^i,^ 



mi + m2 belongs to [0, 2N X — 2], so 



\ 1 „2i7r^(mi+m 2 — m)/N x _ \r cm — mi , at- cjVr+m-ra! 
/ , 31 m 2 ~ x m 2 

H=0 



Similar relations hold for indices n, hence 

T \^ T W fxm—mi , A+m-m^ ( cn-nx , rN y +n— m s 

J-mn — 2__, Jm 1 n 1 VVm 2 n- 2 \°m 2 ' °m 2 J \°n 2 ' °n 2 

m],ni ,m 2 ,ra 2 

Equation 10 specifies that Fourier coefficients are only defined for m, n G [0, N x — 1] X [0, N y — 1] , hence 

rp _ \ rp 1^1,11 

Jmn — / , Jmhil^m! ! 



where 



K 



!W m -m 1 ,n-n 1 if mi < m and iii < n 

W m _ mij jv y+ „- ni if mi < m and ni > 71 

WiVx + m-m! , n-m if ' m l > m and "1 < " 

W Nx ■fm-mi , iVy+n — ?ii ^ m 1 > m an d 711 > 72 



(A.l) 



(A.2) 



(A.3) 



(A.4) 



(A.5) 



(A.6) 



Appendix B: Mask convolution kernel for polarization only 

Polarization maps arc represented in direct space by Stokes parameters Q and U, in angular space by E and B. These parameters are related 

by 
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Qfiv = E [cos(2#£ ni )E mini -sin(2^ ini )B mini ]e 2 ^"^/ JV -+"™i/ Jv v) 

c/ M „ = E W^)^+«(C«)^]« fc<m/,,,MW , 

The polarization pseudo-Fourier coefficients are given by 



N x N y 



V E E E W m 2 - ZCmJ^mm + ^(2^ " Cm,^^] 



mi ,ni m-2 ,n.2 Mi^ 



N X Ny 

>(e 2i7r(pm 1 /N x +mi/N y ) e 2i7v(v,rn.2/N :c +i'n2/Ny) e -2in(fim/N x +vn,/Ny) 

-J— E I- sm(2<C„)Q p „ + cos(2Cr„)^l I^ e -^W".+.'»/",), 

xe 2iir(Mmi /W^+i/ni / N y ) e 2iTT(iim 2 /IV X +vn 2 /N y ) e ~2iT[(p,m / N x + vn/N y ) 

We now have to compute the two summations 

l x — ^^ C Os(20' Jl ' — 2</> M " ^ e 2iTr(pmx/N 3: + vn 1 /Ny) e 2iTr(p.m2/Nx+v7i2/Ny) e -2iTr(pm/Nx+vn/Ny) ^ 

j 2 — ^^ s in(20' Jy — 2</> MI/ } e 2m(ym 1 /Nx+L'ni/Ny) e 2iTr(ij,m2/N x +vn2/Ny) e -2iTv(iim/N x +vn/Ny) _ 



(B.l) 
(B.2) 

(B.3) 
(B.4) 

(B.5) 
(B.6) 

(B.7) 
(B.8) 



Because </> m n is the angle between fe mn and r M „ and <f>m ini is the angle between fc mini and r M „, we have mn — <t>rnxni = ^Ui" ni with 

COs((^' m J 1 1 ) — fcmini ' f^mn • 

As a consequence, the sine and cosine do not depend on pixel indices, fi and u, so we can use the orthogonality relation used for temperature, 
i.e. V mi + wi2 belongs to [0, 2AT — 2] : 

N-l 

E^+^-n)/*. _ J\r rm — ?ni , i\r riV+m — mi 

to finally get 

£ m „ = E ^m',^ 1 bos(2Ci"i)B mi „ 1 +sin(2CJ l ni ) s m 1 n 1 ], (B.9) 
B m „ = £ jr^ [- s in(2Cr i )S mi n 1 +cos(2Cr i )Sm I « I ]. (B.10) 

mi ,ni 

From these last set of results, we can compute the pseudo-power spectra, keeping in mind that 

/ p n p* \ — /-'EE c r 

\ m rn'ri' / ^mn v m : m f u n,n' } 

/ R 71 R* \ (~lBB c r 

\ m-^m'n'/ ^mn u m,m f u n,n' 



and 



We define the estimated pseudo-spectra as 



and 



By using those definitions, we can easily show that 



{ E m B m'n>) ~ { B m E m'n l ) ~ C mn Sm,m> &n,n' 



mn ^mn 



Re -Bmn B * _ — Re B mn E, 



pEE 

mn 


= E 


\ K «,ni |2 


COS 2 (2^™r i )C m fn 1 




m\ ,n\ 






pBB 

mn 


= E 


\ K n,m |2 


)in 2 (20™r i )C m f„ 1 - 




m\ ,n\ 






pEB 


= E 


|7^w,ni |2 
1 v m,jTti 1 


-irin(4#&»i)C£f t 




m\ ,n,\ 
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In a matrix formulation, this reads 



pEE 

mn 
pBB 

mn 
pEB 

mn . 



E 



E 



■ iW m'/j , m ]_ n i 

7i ,rB b ,e e 

MEB,EE 
mn,mi n\ 

nrdiag 
lvl mn,m\ n\ 
fOff 

mn, mi n-i 



ME E . B B 
mn,m\n\ 
71 r B B . B B 
^mn,m\n\ 
nfEB^BB 



EE,EB 
mn,m\ n\ 
BB,EB 
mn,m\ n± 
■EB,EB 
■mn , 777.x n l 




M 



- 2 1V1 ran ,m\n\ 



M °ff 

M off 

lvl mn,m\n-i 
1^ jufcross 
2 mn,m\n 



jurcross \ 
mn,m±n± 

fi/fcross 

mn,m\ n\ 
M diag _ M off I 

lvl mn,m\ n\ lvl mn,m\ni / 



pEE 
pBB 

pEB 

mi ni 



(B.ll) 



(B.12) 



where 



M^% ini = cos 2 (20-r i )l^',^ 1 | 2 , (B.13) 
M°Jlm ini = sin 2 (2^i"i)|AT»^ 1 | 2 , (B.14) 
= sinC^" 1 )!^^! 2 - (B.15) 

When the above mixing matrices are averaged over the two azimuthal (or polar for flat sky) angles, it can be shown that f f M m "°" mini d9d6\ = 
0. However, before such an averaging, it is not a priori zero. 

Appendix C: Mask convolution kernel for temperature polarization cross-correlation 

For the cross-correlation of the temperature with CMB maps, we remind first that 

{T m nE m i n i) = {EmnT^i n i) = C m ^5 mm i8 nn i 

and 

(T mn B m , n ,) = (B mn T m , n ,) = C mn S m m t 8 n n i . 
The estimated cross-pseudo-spectrum are defined as 

Pmn = 2 \_TmnEmn + EmnT mn ^ = Re |^T mn E mn j = Re ^E mn T mn J 

and 

f'mn = ^ \^ mn ^mn + An»^ ml i] = R e ^mn^mn] = -^- e [-^"ii^mn] ■ 
From this and from the above defined pseudo-Fourier coefficients, we show that 



pTE \ . . / M 3aid M^f" \ / P TE 

r mn I \ 1 | lvl mn,m\n\ lvl mn,mini 1 / 1 m\n\ 



pTB I 1 „lf//» M gaid I I pTl 

m l n l V mn,m\n\ J "mn,rn.i»li y \ mini 



(ci) 



with 

M^i d mini = cos(2^r i )l^ 1 | 2 , (C2) 
M m f n ° mini = sin(2^r i )l^™'^ 1 | 2 - (C3) 

As for the cross blocks in the polarization case, the azimuthal average of Mmn°,min\ vanishes, i.e. J J M£/ n ° mini d9d9i = 0. However, before 
this averaging, it is not a priori zero. 
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